/* José Filipe M Santos - 7277502 */

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <gmp.h>
#include <math.h>
#include <time.h>

#define DIGITOS 10000000

int i;

/* execucao do algoritmo de borwein com convergencia quadratica */
void borwein(){

	i = 0; /* quantidade de iteracoes */

	unsigned long int um = 1, dois = 2; /* ctes utilizadas na formula */

	mpf_t a0, b0, p0, sqrta0, intermedio1, intermedio2, pi, diff, zero; 	/* a0, b0, p0: variaveis da formula 				*/
										/* srtqa0: raiz de a0     					*/
										/* intermedio[i]: valores intermediarios utilizados nos calculos*/
										/* pi, diff, zero: valores utilizados para verificar cvergencia	*/

	FILE *ifp = fopen("pi.txt", "w+");

	mpf_set_default_prec(DIGITOS * sizeof(int));				/* precisao de 10 mi de digitos */
	mpf_init(a0);								/* inicializando floats		*/
	mpf_init(sqrta0);
	mpf_init(b0);
	mpf_init(p0);
	mpf_init(pi);
	mpf_init(diff);
	mpf_init(intermedio1);
	mpf_init(intermedio2);
	mpf_init(zero);

	mpf_sqrt_ui(a0, dois);							/* fixando valores iniciais (a0 = sqrt(2) e p0 = 2 + sqrt(2))*/
	mpf_sqrt_ui(p0, dois);
	mpf_add_ui(p0, p0, dois);

	do{

		mpf_set(pi, p0);						/* pi = p0 atual*/
		mpf_sqrt(sqrta0, a0);

		mpf_add_ui(intermedio1, b0, um); 				/* int1 = (1 + b0) */
		mpf_mul(intermedio1, intermedio1, sqrta0); 			/* int1 = (sqrt(a0)) * (1 + b0) */
		mpf_add(intermedio2, a0, b0); 					/* int2 = a0 + b0 */
		mpf_div(b0, intermedio1, intermedio2); 				/* b0 = (sqrt(a0) * (1 + b0)) / (a0 + b0) */

										/* b[n+1] = (sqrt(a[n]) * (1 + b[n])) / (a[n] + bn[n]) */

		mpf_add_ui(intermedio1, a0, um);
		mpf_mul_ui(intermedio2, sqrta0, dois);
		mpf_div(a0, intermedio1, intermedio2);

										/* a[n+1] = a[n] + 1 / 2 x sqrt(a[n]) */

		mpf_add_ui(intermedio1, a0, 1);
		mpf_add_ui(intermedio2, b0, 1);
		mpf_div(intermedio2, intermedio1, intermedio2);
		mpf_mul(intermedio1, b0, p0);
		mpf_mul(p0, intermedio1, intermedio2);

										/* p[n+1] = (p[n] * b[n+1])*((1 + a[n+1])/(1 + b[n+1])) */

		mpf_sub(diff, p0, pi);						/* tirar a diferenca, se for 0, está pronto */
		++i;
		mpf_out_str (ifp, 10, DIGITOS, p0);				/* inserir pi a cada iteracao no arquivo */
		fprintf(ifp, " \n");

	}while(mpf_cmp(zero, diff) != 0);					/* iterar ate diff == 0 */

	fprintf(ifp, " iteracoes: %d\n", i);					/* inserir numero total de iteracoes */
	fclose(ifp);

	return;

}

int main(){

	time_t inicio, fim;

	time(&inicio);
	borwein();
	time(&fim);
	printf("RELATORIO DE EXECUCAO:");
	printf("%d ITERACOES EM %lfs\n", i, difftime(fim, inicio));

	return 0;
}
